A synthesis of dryland restoration techniques.

Purpose

To quantitatively examine the efficacy of vegetation restoration in drylands globally.

Questions

  1. What is the global extent of research that directly examined restoration of drylands?
  2. What were the common measures?
  3. Is the restoration of vegetation a common and primary focus?
  4. How frequently does the restoration measure outcomes beyond the focal species?
  5. What were the primary restoration goals as reported by primary authors?
  6. How much variation was there in the techniques tested and how long were experiments monitored and tested?
  7. How relatively effective were the techniques?

Step 2. Sort

A summary of sort process using PRISMA.

library(PRISMAstatement)
prisma(found = 1504,
       found_other = 5,
       no_dupes = 1039, 
       screened = 1039, 
       screen_exclusions = 861, 
       full_text = 178,
       full_text_exclusions = 101, 
       qualitative = 77, 
       quantitative = 66,
       width = 800, height = 800)

Step 3. Synthesize

Check data and calculate necessary measures.

#all data
#data <- read_csv("data/data.csv")
#data <- data %>%
  #mutate(lrr = log(mean.t/mean.c), rii = ((mean.t-mean.c)/(mean.t + mean.c)), var.es = ((sd.t^2/n.t*mean.t^2) + (sd.c^2/n.c*mean.c^2)))

#other effect size estimates
#library(compute.es)
#data <- data %>%
  #mutate(d=mes(mean.t, mean.c, sd.t, sd.c, n.t, n.c, level = 95, #cer = 0.2, dig = 2, , id = ID, data = data))
#check metafor

#data from ag and grazing studies that examined restoration in drylands
data <- read_csv("data/data.csv")
mydata<- data %>% 
 filter(disturbance %in% c("agriculture","grazing")) %>% 
filter(!notes %in% "couldnt extract data")

#write.csv(mydata, file="mydata.csv")

mydata <- mydata %>%
  mutate(lrr = log(mean.t/mean.c), rii = ((mean.t-mean.c)/(mean.t + mean.c)), var.es = ((sd.t^2/(n.t*mean.t^2)) + (sd.c^2/(n.c*mean.c^2)))) #use only lrr now

#mydata<- mydata %>% 
 # mutate(aridity.range = cut(aridity.index, breaks=c(0,5,10,20,25,69.92), labels=c("hyperarid","arid", #"semiarid", "moderately arid", "slightly humid"))) #categorization of aridity.index values, according to de Martone

Step 4. Summarize

Explore summary level data of all data. Explore aggregation levels that support the most reasonable data structure and minimize non-independence issues.

#evidence map####
require(maps)
world<-map_data("world")
map<-ggplot() + geom_polygon(data=world, fill="gray50", aes(x=long, y=lat, group=group))
#map + geom_point(data=paperdata, aes(x=long, y=lat)) +
  #labs(x = "longitude", y = "latitude") #render a literal map, i.e. evidence map, of where we study the niche in deserts globally

#add in levels and color code points on map####
#map of 178 articles
map + geom_point(data=data, aes(x=long, y=lat, color = paradigm)) + 
  scale_color_brewer(palette = "Paired") +
  labs(x = "longitude", y = "latitude", color = "")

#map of selected articles, agriculture and grazing disturbances
map + geom_point(data=mydata, aes(x=long, y=lat, color = paradigm)) + 
  scale_color_brewer(palette = "Paired") +
  labs(x = "longitude", y = "latitude", color = "")

#aggregation####
se <- function(x){
  sd(x)/sqrt(length(x))
}

data.simple <- mydata %>%
  group_by(study.ID, paradigm, technique, measure.success) %>%
  summarise(n = n(), mean.lrr = mean(lrr), mean.rii = mean(rii), mean.var = mean(var.es))

main.data <- mydata %>%
  group_by(study.ID, paradigm, intervention, outcome) %>%
  summarise(n = n(), mean.lrr = mean(lrr), mean.rii = mean(rii), mean.var = mean(var.es))


#EDA data####
simple.data <- mydata %>% group_by(study.ID, paradigm, technique, measure.success) %>% summarise(mean.rii = mean(rii), error = se(rii))
simple.data <- na.omit(simple.data)

parad.data <- mydata %>% group_by(study.ID, paradigm) %>% summarise(mean.rii = mean(rii), error = se(rii))
parad.data <- na.omit(parad.data)

tech.data <- mydata %>% group_by(study.ID, technique) %>% summarise(mean.rii = mean(rii), error = se(rii))
tech.data <- na.omit(tech.data)

success.data <- mydata %>% group_by(study.ID, measure.success) %>% summarise(mean.rii = mean(rii), error = se(rii))
success.data <- na.omit(success.data)


#active
active <- mydata %>%
  filter(paradigm == "active")

#viz for aggregation####
disturbance.data <- mydata %>% 
  group_by(study.ID,disturbance, paradigm) %>% 
  count()

disturbance.data2 <- disturbance.data %>% 
  group_by(disturbance,paradigm) %>% 
  count()

ggplot(na.omit(disturbance.data2), aes(disturbance,nn, fill=paradigm))+ 
  geom_bar(stat = "identity") + 
  coord_flip(ylim=0:44) + 
  scale_fill_brewer(palette = "Blues")

intervention.data <- active %>% 
  group_by(study.ID,intervention, paradigm) %>% 
  count()

intervention.data2 <- intervention.data %>% 
  group_by(intervention,paradigm) %>% 
  count()

#ggplot(na.omit(intervention.data2), aes(intervention,nn, fill=paradigm)) +
  #geom_bar(stat = "identity") + 
  #coord_flip(ylim=0:44) + 
  #scale_fill_brewer(palette = "Blues")

technique.data <- mydata %>% 
  group_by(study.ID,technique, paradigm) %>% 
  count()

technique.data2 <- technique.data %>% 
  group_by(technique,paradigm) %>% 
  count()

technique.data3 <- technique.data2 %>% 
  group_by(technique) %>% 
  count()

aridity<- mydata %>% 
  group_by(continent) %>% summarize(mean.aridity=mean(aridity.index))
aridity  
## # A tibble: 6 x 2
##   continent     mean.aridity
##   <chr>                <dbl>
## 1 Africa                9.04
## 2 Asia                 20.7 
## 3 Europe               26.8 
## 4 North America        19.6 
## 5 Oceania              16.4 
## 6 South America        22.9
#table(mydata$aridity.index)


ggplot(na.omit(data.simple), aes(technique, n, fill = paradigm)) + 
  geom_bar(stat = "identity") + 
  coord_flip() + 
  scale_fill_brewer(palette = "Paired")

ggplot(na.omit(data.simple), aes(measure.success, n, fill = paradigm)) + 
  geom_bar(stat = "identity") + 
  coord_flip() + 
  scale_fill_brewer(palette = "Paired")

#better
ggplot(main.data, aes(intervention, n, fill = paradigm)) + 
  geom_bar(stat = "identity") + 
  coord_flip() + 
  scale_fill_brewer(palette = "Paired") +
  labs(fill = "")

ggplot(main.data, aes(outcome, n, fill = paradigm)) + 
  geom_bar(stat = "identity") + 
  coord_flip() + 
  scale_fill_brewer(palette = "Paired") +
  labs(fill = "")

Step 5a. EDA

Exploratory data analyses to understand data and QA/QC using Rii.

Step 5b. Models

Meta and conventional statistical models to explore relative efficacy.

Step 5. Synthesis stats

#p-value meta
#library(metap)
#all data for metas but cleaned for na's
#all data for meta cleaned
mdata.all <- mydata %>%
  filter(!is.na(lrr)) %>%
  filter(!is.na(var.es)) %>%
  filter(!is.na(n.t)) %>%
  filter(!is.na(p)) %>%
  filter(!is.na(intervention)) %>%
  filter(is.finite(lrr)) %>%
  filter(!is.na(exp.length)) %>%
  filter(!is.na(MAP)) %>%
  filter(!is.na(aridity.index))
 
#after filtering NAs we apply models to 40 articles

summary(mdata.all$var.es)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.00000  0.00000  0.00076  0.15430  0.01441 39.06266
#plot(density(mdata.all$var.es))
top_n(mdata.all, -10, var.es)
## # A tibble: 10 x 50
##    study.ID    ID publication.year data.year exp.year exp.length
##       <dbl> <dbl>            <dbl>     <dbl>    <dbl>      <dbl>
##  1      111  111.             2016      2003     2003        132
##  2      111  111.             2016      2003     2003        132
##  3      111  111.             2016      2003     2003        132
##  4      111  111.             2016      2003     2003        132
##  5      111  111.             2016      2003     2003        132
##  6      111  111.             2016      2003     2003        132
##  7      111  111.             2016      2003     2003        132
##  8      111  111.             2016      2003     2003        132
##  9      111  111.             2016      2003     2003        132
## 10      111  111.             2016      2003     2003        132
## # ... with 44 more variables: disturbance <chr>, focus <chr>,
## #   technique <chr>, intervention <chr>, paradigm <chr>, hypothesis <chr>,
## #   pathway <chr>, plant.species <chr>, target.plant <chr>,
## #   measure.success <chr>, outcome <chr>, measured.factor <chr>,
## #   treatment <chr>, control <chr>, unit <chr>, Nsites <dbl>, n.t <dbl>,
## #   n.c <dbl>, ntotalsamples <dbl>, mean.t <dbl>, mean.c <dbl>,
## #   sd.t <dbl>, sd.c <dbl>, se.t <dbl>, se.c <dbl>, p <dbl>, df <dbl>,
## #   measure.dispersion <chr>, lat <dbl>, long <dbl>, continent <chr>,
## #   country <chr>, ecosystem <chr>, elevation <dbl>, MAP <dbl>, MAT <dbl>,
## #   aridity.index <dbl>, `potential.evaporation (mm)` <dbl>,
## #   grazing <chr>, soil <chr>, notes <chr>, lrr <dbl>, rii <dbl>,
## #   var.es <dbl>
top_n(mdata.all, 5, var.es)
## # A tibble: 5 x 50
##   study.ID    ID publication.year data.year exp.year exp.length disturbance
##      <dbl> <dbl>            <dbl>     <dbl>    <dbl>      <dbl> <chr>      
## 1       11  11.7             2015      2007     1956        600 agriculture
## 2       11  11.7             2015      2007     1956        600 agriculture
## 3       11  11.7             2015      2007     1956        600 agriculture
## 4       11  11.7             2015      2007     1956        600 agriculture
## 5      147 147.              2012      2005     2005         36 agriculture
## # ... with 43 more variables: focus <chr>, technique <chr>,
## #   intervention <chr>, paradigm <chr>, hypothesis <chr>, pathway <chr>,
## #   plant.species <chr>, target.plant <chr>, measure.success <chr>,
## #   outcome <chr>, measured.factor <chr>, treatment <chr>, control <chr>,
## #   unit <chr>, Nsites <dbl>, n.t <dbl>, n.c <dbl>, ntotalsamples <dbl>,
## #   mean.t <dbl>, mean.c <dbl>, sd.t <dbl>, sd.c <dbl>, se.t <dbl>,
## #   se.c <dbl>, p <dbl>, df <dbl>, measure.dispersion <chr>, lat <dbl>,
## #   long <dbl>, continent <chr>, country <chr>, ecosystem <chr>,
## #   elevation <dbl>, MAP <dbl>, MAT <dbl>, aridity.index <dbl>,
## #   `potential.evaporation (mm)` <dbl>, grazing <chr>, soil <chr>,
## #   notes <chr>, lrr <dbl>, rii <dbl>, var.es <dbl>
max(mdata.all$var.es)/min(mdata.all$var.es)
## [1] 4.800933e+14
table(mdata.all$study.ID)
## 
##  11  13  51  64  66  68  69  77  87  88  95  98 104 109 111 112 115 119 
##  84  26  24  45  24  12  10   3  48 117  15   2 128   6 249   7  21  24 
## 121 135 142 147 152 154 161 167 170 179 206 207 210 223 231 239 247 251 
## 108  30  12  32   7  12   4  18   3  12  63  15  27  15   7   4 198   8 
## 255 256 263 264 
##   5   6  20   9
#active only data for meta cleaned up
mdata <- mydata %>%
  filter(paradigm == "active") %>%
  filter(!is.na(lrr)) %>%
  filter(!is.na(var.es)) %>%
  filter(!is.na(n.t)) %>%
  filter(!is.na(p)) %>%
  filter(!is.na(intervention)) %>%
  filter(is.finite(lrr)) %>%
  filter(!is.na(exp.length))

#passive only data for meta cleaned up
mdatap <- mydata %>%
  filter(paradigm == "passive") %>%
  filter(!is.na(lrr)) %>%
  filter(!is.na(var.es)) %>%
  filter(!is.na(n.t)) %>%
  filter(!is.na(p)) %>%
  filter(!is.na(intervention)) %>%
  filter(is.finite(lrr)) %>%
  filter(!is.na(exp.length))

#aggregated active data for metas var estimated with central tendency
#note - could also bootstrap mean variance here instead of arithematic mean
simple.mdata <- mdata %>%
  group_by(intervention) %>%
  summarise(lrr = mean(lrr), var.es = mean(var.es), n = mean(n.t))

simple.mdata.2 <- mdata %>%
  group_by(intervention, outcome) %>%
  summarise(lrr = mean(lrr), var.es = mean(var.es), n = mean(n.t))  

simple.mdata.var <- mdata %>%
  group_by(intervention) %>%
  summarise(lrr = mean(lrr), var.es = se(var.es), n = mean(n.t))

simple.mdata2.var <- mdata %>%
  group_by(intervention, outcome) %>%
  summarise(lrr = mean(lrr), var.es = se(var.es), n = mean(n.t))


#aggregated passive data for metas var estimated with central tendency
simple.mdatap <- mdatap %>%
  group_by(intervention) %>%
  summarise(lrr = mean(lrr), var.es = mean(var.es), n = mean(n.t))

simple.mdatap.2 <- mdatap %>%
  group_by(intervention, outcome) %>%
  summarise(lrr = mean(lrr), var.es = mean(var.es), n = mean(n.t))  

simple.mdatap.var <- mdatap %>%
  group_by(intervention) %>%
  summarise(lrr = mean(lrr), var.es = se(var.es), n = mean(n.t))

simple.mdatap2.var <- mdatap %>%
  group_by(intervention, outcome) %>%
  summarise(lrr = mean(lrr), var.es = se(var.es), n = mean(n.t))


#metas with p-values####
#schweder(mdata$p)
#sumz(p, data = mdata)
#mdata %>%
  #split(.$intervention) %>%
  #purrr::map(~sumz(.$p)) 
#sumlog(mdata$p)

#metas with meta package on effect size measures####
library(meta)
#all data non-aggregated
#m <- metagen(lrr, var.es, studlab = ID, byvar = intervention, data = mdata)
#summary(m)
#funnel(m)
#metabias(m)
#forest(m, sortvar = intervention)

#t-tests if different from 0

tmu <- function(x){t.test(x, mu = 0, paired = FALSE, var.equal=FALSE, conf.level = 0.95)}

mdata.all %>%
  split(.$paradigm) %>%
  purrr::map(~tmu(.$lrr)) #note this uses arithmetic means not estimated means from random effect models
## $active
## 
##  One Sample t-test
## 
## data:  x
## t = 7.6083, df = 1101, p-value = 5.943e-14
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.2069479 0.3507825
## sample estimates:
## mean of x 
## 0.2788652 
## 
## 
## $passive
## 
##  One Sample t-test
## 
## data:  x
## t = -7.5438, df = 357, p-value = 3.824e-13
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.4408271 -0.2585133
## sample estimates:
##  mean of x 
## -0.3496702
#Propose we only use random effect models because of the diversity of studies, differences in the methods and samples that may introduce heterogeneity
m1 <- metagen(lrr, var.es, studlab = ID, byvar = paradigm, comb.fixed=FALSE, data = mdata.all)
summary(m1) #active is positive and passive is negative so should not mix
## Number of studies combined: k = 1460
## 
##                                       95%-CI     z  p-value
## Random effects model 0.0766 [0.0654; 0.0879] 13.38 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 184527270028.74 [184527270027.80; 184527270029.68]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                              Q d.f. p-value
##  49679407227636228414242816.00 1459       0
## 
## Results for subgroups (random effects model):
##                       k                     95%-CI
## paradigm = active  1102  0.2184 [ 0.2055;  0.2314]
## paradigm = passive  358 -0.3413 [-0.3753; -0.3073]
##                                                Q  tau^2    I^2
## paradigm = active  49671693556322550581035008.00 0.0450 100.0%
## paradigm = passive     4152232320954931871744.00 0.1047 100.0%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   911.23    1 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
mb1<- metabias(m1)
mb1
#should do a t-test right now of paradigm is diff from 0
metareg(m1,aridity.index+exp.length)
## 
## Mixed-Effects Model (k = 1460; tau^2 estimator: DL)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0440 (SE = 0.0167)
## tau (square root of estimated tau^2 value):             0.2098
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   33284207035718112903168.00
## R^2 (amount of heterogeneity accounted for):            2.22%
## 
## Test for Residual Heterogeneity:
## QE(df = 1457) = 48495089651041292387352576.0000, p-val < .0001
## 
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 3535.8262, p-val < .0001
## 
## Model Results:
## 
##                estimate      se      zval    pval    ci.lb    ci.ub 
## intrcpt          0.8886  0.0150   59.3593  <.0001   0.8593   0.9179  *** 
## aridity.index   -0.0339  0.0007  -51.0096  <.0001  -0.0352  -0.0326  *** 
## exp.length      -0.0007  0.0000  -26.2679  <.0001  -0.0008  -0.0006  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#interventions
m2 <- metagen(lrr, var.es, studlab = ID, byvar = intervention, comb.fixed=FALSE, subset = paradigm == "active", data = mdata.all)
summary(m2)
## Number of studies combined: k = 1102
## 
##                                       95%-CI     z  p-value
## Random effects model 0.2184 [0.2055; 0.2314] 33.02 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 212403086959.62 [212403086958.53; 212403086960.71]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                              Q d.f. p-value
##  49671693556322550581035008.00 1101       0
## 
## Results for subgroups (random effects model):
##                                 k                  95%-CI
## intervention = vegetation     779 0.1845 [0.1694; 0.1996]
## intervention = soil           248 0.3128 [0.2990; 0.3265]
## intervention = water addition  75 0.6409 [0.5539; 0.7279]
##                                                           Q  tau^2    I^2
## intervention = vegetation     48393423300974332080029696.00 0.0443 100.0%
## intervention = soil              97513761686148203151360.00 0.0111 100.0%
## intervention = water addition                      31132.18 0.1047  99.8%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   226.58    2 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
#metabias(m)
metareg(m2, ~ aridity.index + exp.length)
## 
## Mixed-Effects Model (k = 1102; tau^2 estimator: DL)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0440 (SE = 0.0167)
## tau (square root of estimated tau^2 value):             0.2098
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   44122812844877792935936.00
## R^2 (amount of heterogeneity accounted for):            2.22%
## 
## Test for Residual Heterogeneity:
## QE(df = 1099) = 48490971316520689160159232.0000, p-val < .0001
## 
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 1836.2421, p-val < .0001
## 
## Model Results:
## 
##                estimate      se      zval    pval    ci.lb    ci.ub 
## intrcpt          0.5607  0.0193   29.1251  <.0001   0.5230   0.5985  *** 
## aridity.index   -0.0237  0.0008  -29.5763  <.0001  -0.0253  -0.0221  *** 
## exp.length       0.0009  0.0000   19.0317  <.0001   0.0008   0.0010  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m3 <- metagen(lrr, var.es, studlab = ID, byvar = intervention, subset = paradigm == "passive", comb.fixed=FALSE, data = mdata.all)
summary(m3)
## Number of studies combined: k = 358
## 
##                                          95%-CI      z  p-value
## Random effects model -0.3413 [-0.3753; -0.3073] -19.70 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.1047; H = 3410410951.75 [3410410950.14; 3410410953.36]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                          Q d.f. p-value
##  4152232320954931871744.00  357       0
## 
## Results for subgroups (random effects model):
##                                    k                     95%-CI
## intervention = vegetation        125  0.2654 [ 0.2067;  0.3241]
## intervention = grazing exclusion  29  0.1351 [ 0.0270;  0.2431]
## intervention = soil              204 -0.7583 [-0.8196; -0.6970]
##                                                          Q  tau^2    I^2
## intervention = vegetation        4152209524073903423488.00 0.1047 100.0%
## intervention = grazing exclusion              238316232.18 0.0881 100.0%
## intervention = soil                   14453104123321616.00 0.1990 100.0%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   595.91    2 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
#metabias(m)
metareg(m3, ~ aridity.index + exp.length)
## 
## Mixed-Effects Model (k = 358; tau^2 estimator: DL)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.1047 (SE = 0.0876)
## tau (square root of estimated tau^2 value):             0.3236
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   11696393769174456320.00
## R^2 (amount of heterogeneity accounted for):            0.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 355) = 4152219788056932122624.0000, p-val < .0001
## 
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 815.4196, p-val < .0001
## 
## Model Results:
## 
##                estimate      se      zval    pval    ci.lb    ci.ub 
## intrcpt          0.3781  0.0565    6.6935  <.0001   0.2674   0.4888  *** 
## aridity.index    0.0015  0.0029    0.5063  0.6127  -0.0042   0.0072      
## exp.length      -0.0020  0.0001  -23.0381  <.0001  -0.0021  -0.0018  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#trying models by continents
#m4 <- metagen(lrr, var.es, studlab = ID, byvar = continent, subset = paradigm == "active", comb.fixed=FALSE, data = #mdata.all)
#summary(m4)

#m5 <- metagen(lrr, var.es, studlab = ID, byvar = continent, subset = paradigm == "passive", comb.fixed=FALSE, data = #mdata.all)
#summary(m5)


#use random effects mean and var estimate to plot

#do t-tests here too

#outcomes
m6 <- metagen(lrr, var.es, studlab = ID, byvar = outcome, subset = paradigm == "active", comb.fixed=FALSE, data = mdata.all)
summary(m6)
## Number of studies combined: k = 1102
## 
##                                       95%-CI     z  p-value
## Random effects model 0.2184 [0.2055; 0.2314] 33.02 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.0450; H = 212403086959.62 [212403086958.53; 212403086960.71]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                              Q d.f. p-value
##  49671693556322550581035008.00 1101       0
## 
## Results for subgroups (random effects model):
##                     k                     95%-CI
## outcome = soil    249  0.2204 [ 0.1558;  0.2849]
## outcome = plants  305  0.5071 [ 0.4936;  0.5206]
## outcome = animals  24 -0.1152 [-0.1155; -0.1148]
## outcome = habitat 524  0.0621 [ 0.0437;  0.0804]
##                                               Q   tau^2    I^2
## outcome = soil                35077220764051.67  0.2656 100.0%
## outcome = plants     97513760543782884868096.00  0.0111 100.0%
## outcome = animals                     541696.64 <0.0001 100.0%
## outcome = habitat 48393423303339003634253824.00  0.0443 100.0%
## 
## Test for subgroup differences (random effects model):
##                        Q d.f. p-value
## Between groups   8647.81    3       0
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
#metabias(m)
metareg(m6, ~ aridity.index + exp.length)
## 
## Mixed-Effects Model (k = 1102; tau^2 estimator: DL)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0440 (SE = 0.0167)
## tau (square root of estimated tau^2 value):             0.2098
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   44122812844877792935936.00
## R^2 (amount of heterogeneity accounted for):            2.22%
## 
## Test for Residual Heterogeneity:
## QE(df = 1099) = 48490971316520689160159232.0000, p-val < .0001
## 
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 1836.2421, p-val < .0001
## 
## Model Results:
## 
##                estimate      se      zval    pval    ci.lb    ci.ub 
## intrcpt          0.5607  0.0193   29.1251  <.0001   0.5230   0.5985  *** 
## aridity.index   -0.0237  0.0008  -29.5763  <.0001  -0.0253  -0.0221  *** 
## exp.length       0.0009  0.0000   19.0317  <.0001   0.0008   0.0010  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m7 <- metagen(lrr, var.es, studlab = ID, byvar = outcome, subset = paradigm == "passive", comb.fixed=FALSE, data = mdata.all)
summary(m7)
## Number of studies combined: k = 358
## 
##                                          95%-CI      z  p-value
## Random effects model -0.3413 [-0.3753; -0.3073] -19.70 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 0.1047; H = 3410410951.75 [3410410950.14; 3410410953.36]; I^2 = 100.0% [100.0%; 100.0%]
## 
## Test of heterogeneity:
##                          Q d.f. p-value
##  4152232320954931871744.00  357       0
## 
## Results for subgroups (random effects model):
##                     k                     95%-CI                         Q
## outcome = habitat 104  0.1605 [ 0.0964;  0.2246] 4152172019980636258304.00
## outcome = plants   50  0.4438 [ 0.0345;  0.8532]      33314950066906688.00
## outcome = soil    204 -0.7583 [-0.8196; -0.6970]      14453104123321616.00
##                    tau^2    I^2
## outcome = habitat 0.1047 100.0%
## outcome = plants  2.1620 100.0%
## outcome = soil    0.1990 100.0%
## 
## Test for subgroup differences (random effects model):
##                       Q d.f.  p-value
## Between groups   425.72    2 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
#metabias(m)
metareg(m7, ~ aridity.index + exp.length)
## 
## Mixed-Effects Model (k = 358; tau^2 estimator: DL)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.1047 (SE = 0.0876)
## tau (square root of estimated tau^2 value):             0.3236
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   11696393769174456320.00
## R^2 (amount of heterogeneity accounted for):            0.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 355) = 4152219788056932122624.0000, p-val < .0001
## 
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 815.4196, p-val < .0001
## 
## Model Results:
## 
##                estimate      se      zval    pval    ci.lb    ci.ub 
## intrcpt          0.3781  0.0565    6.6935  <.0001   0.2674   0.4888  *** 
## aridity.index    0.0015  0.0029    0.5063  0.6127  -0.0042   0.0072      
## exp.length      -0.0020  0.0001  -23.0381  <.0001  -0.0021  -0.0018  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#t-tests for outcomes diff from 0 but you can see using the 95% CI what is different

#check metafor for interactions?? in these big models or are we ok??

#brainstorm on how to explore?? techniques


#save but cut all this.
#m.study <- metagen(lrr, var.es, studlab = study.ID, byvar = intervention, data = mdata)
#summary(m.study)
#funnel(m)
#metabias(m)
#forest(m, sortvar = intervention)

#aggregated data
#m1 <- metagen(lrr, var.es, studlab = intervention, data = simple.mdata)
#summary(m1)
#funnel(m1)
#metabias(m1)
#forest(m, sortvar = intervention)

#different var estimate
#m2 <- metagen(lrr, var.es, studlab = intervention, data = simple.mdata.var)
#summary(m2)
#funnel(m2)
#metabias(m2)
#forest(m, sortvar = intervention)


#m3 <- metagen(lrr, var.es, studlab = intervention, byvar = outcome, data = simple.mdata.2)
#summary(m3)
#funnel(m)
#radial(m3)
#metabias(m2)
#forest(m, sortvar = intervention)

#m4 <- metagen(lrr, var.es, studlab = intervention, byvar = outcome, data = simple.mdata2.var)
#summary(m4)
#funnel(m)
#radial(m4)
#metabias(m2)
#forest(m, sortvar = intervention)
library(metafor)

#With metafor it is possible to fit metaregression models (e.g., Berkey et al. 1995; van Houwelingen et al. 2002), that #is, linear models that examine the infuence of one or more moderator variables on the outcomes (Viechtbauer, 2010). 

# I´ve installed the 'devel' version of metafor (https://wviechtb.github.io/metafor/#installation)

#checking and comparing Lrr and sampling variances values obtained with escalc() function and by functions manually defined (lines 153-155)

#metadata<-escalc(measure="ROM",m1i=mean.t,m2i=mean.c,sd1i=sd.t,sd2i=sd.c,n1i=n.t,n2i=n.c, #data=mydata,var.names=c("LRR","LRR_var"),digits=4)


#metadata <- metadata %>%
#  filter(!is.na(LRR)) %>%
#  filter(!is.na(LRR_var)) %>%
#  filter(!is.na(n.t)) %>%
#  filter(!is.na(p)) %>%
#  filter(!is.na(intervention)) %>%
#  filter(is.finite(lrr)) %>%
#  filter(!is.na(exp.length)) %>% 
#  filter(!is.na(aridity.index))

#summary(metadata$LRR_var)
#plot(density(metadata$var.es))
#top_n(metadata, -10, LRR_var)
#top_n(metadata, 5, LRR_var)

#write.csv(metadata, "metadata.csv")

#mdata.veg <- mydata %>%
 # filter(paradigm == "active") %>%
  #filter(intervention=="vegetation") %>% #filtering by intervention
  #filter(!is.na(lrr)) %>%
  #filter(!is.na(var.es)) %>%
#  filter(!is.na(n.t)) %>%
#  filter(!is.na(p)) %>%
#  filter(!is.na(intervention)) %>%
#  filter(is.finite(lrr)) %>%
#  filter(!is.na(exp.length))

mdata.mean<- mdata.all %>% 
  group_by(study.ID, aridity.index, exp.length, paradigm, intervention, outcome) %>% 
  summarize(mean_lrr= mean(lrr),
  mean_var.es= mean(var.es)) #collapsing data to means

#mdata.mean<- mdata.mean %>%
#  filter(!is.infinite(mean_var.es)) #filter out outliers that are infinite


mod.1 <- rma(yi=lrr, vi=var.es, data = mdata.all)
summary(mod.1)
## 
## Random-Effects Model (k = 1460; tau^2 estimator: REML)
## 
##     logLik    deviance         AIC         BIC        AICc 
## -2057.9849   4115.9697   4119.9697   4130.5407   4119.9780   
## 
## tau^2 (estimated amount of total heterogeneity): 0.8959 (SE = 0.0341)
## tau (square root of estimated tau^2 value):      0.9465
## I^2 (total heterogeneity / total variability):   100.00%
## H^2 (total variability / sampling variability):  96407346870.76
## 
## Test for Heterogeneity:
## Q(df = 1459) = 7665324384897.3271, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub 
##   0.0899  0.0252  3.5678  0.0004  0.0405  0.1393  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(0.0899) #estimate
## [1] 1.094065
mod.2 <- rma(lrr, var.es, mods=  ~paradigm+aridity.index+exp.length -1, data = mdata.all)
summary(mod.2)
## 
## Mixed-Effects Model (k = 1460; tau^2 estimator: REML)
## 
##     logLik    deviance         AIC         BIC        AICc 
## -1928.4827   3856.9655   3866.9655   3893.3827   3867.0069   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.7501 (SE = 0.0287)
## tau (square root of estimated tau^2 value):             0.8661
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   80126342989.09
## 
## Test for Residual Heterogeneity:
## QE(df = 1456) = 6538151352462.4238, p-val < .0001
## 
## Test of Moderators (coefficients 1:4):
## QM(df = 4) = 290.7440, p-val < .0001
## 
## Model Results:
## 
##                  estimate      se      zval    pval    ci.lb    ci.ub 
## paradigmactive     0.9059  0.0601   15.0662  <.0001   0.7880   1.0237  *** 
## paradigmpassive    0.5340  0.0924    5.7816  <.0001   0.3529   0.7150  *** 
## aridity.index     -0.0337  0.0027  -12.4170  <.0001  -0.0391  -0.0284  *** 
## exp.length        -0.0003  0.0001   -2.2753  0.0229  -0.0005  -0.0000    * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Evaluating the influence of moderators on the heterogeneity

#interventions
mod.3 <- rma(lrr, var.es, slab= ID, mods=  ~intervention+aridity.index+exp.length -1, data = mdata.all, subset = paradigm == "active")
summary(mod.3)
## 
## Mixed-Effects Model (k = 1102; tau^2 estimator: REML)
## 
##     logLik    deviance         AIC         BIC        AICc 
## -1470.7522   2941.5045   2953.5045   2983.5065   2953.5816   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.7689 (SE = 0.0339)
## tau (square root of estimated tau^2 value):             0.8769
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   108731387818.35
## 
## Test for Residual Heterogeneity:
## QE(df = 1097) = 6497632603139.7705, p-val < .0001
## 
## Test of Moderators (coefficients 1:5):
## QM(df = 5) = 237.3672, p-val < .0001
## 
## Model Results:
## 
##                             estimate      se     zval    pval    ci.lb 
## interventionsoil              0.7767  0.0959   8.0997  <.0001   0.5887 
## interventionvegetation        0.4676  0.0869   5.3814  <.0001   0.2973 
## interventionwater addition    0.9680  0.1197   8.0847  <.0001   0.7333 
## aridity.index                -0.0258  0.0034  -7.5588  <.0001  -0.0324 
## exp.length                    0.0013  0.0002   5.9590  <.0001   0.0009 
##                               ci.ub 
## interventionsoil             0.9646  *** 
## interventionvegetation       0.6379  *** 
## interventionwater addition   1.2027  *** 
## aridity.index               -0.0191  *** 
## exp.length                   0.0017  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#forest(mod.3, slab= "study.ID")

mod.4 <- rma(lrr, var.es, slab= ID, mods=  ~intervention+aridity.index+exp.length -1, data = mdata.all, subset = paradigm == "passive")
summary(mod.4) #aridity.index and exp.length non significants
## 
## Mixed-Effects Model (k = 358; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc 
## -375.9568   751.9137   763.9137   787.1125   764.1565   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.4459 (SE = 0.0350)
## tau (square root of estimated tau^2 value):             0.6677
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   498110093.57
## 
## Test for Residual Heterogeneity:
## QE(df = 353) = 37647543864.7879, p-val < .0001
## 
## Test of Moderators (coefficients 1:5):
## QM(df = 5) = 275.3915, p-val < .0001
## 
## Model Results:
## 
##                                estimate      se     zval    pval    ci.lb 
## interventiongrazing exclusion    0.2630  0.1818   1.4465  0.1480  -0.0934 
## interventionsoil                -0.1622  0.3593  -0.4515  0.6516  -0.8664 
## interventionvegetation           0.3119  0.1185   2.6329  0.0085   0.0797 
## aridity.index                    0.0010  0.0064   0.1495  0.8812  -0.0116 
## exp.length                      -0.0010  0.0007  -1.5800  0.1141  -0.0023 
##                                 ci.ub 
## interventiongrazing exclusion  0.6194     
## interventionsoil               0.5420     
## interventionvegetation         0.5441  ** 
## aridity.index                  0.0135     
## exp.length                     0.0003     
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.4b <- rma(lrr, var.es, slab= ID, mods=  ~intervention -1, data = mdata.all, subset = paradigm == "passive")
summary(mod.4b)
## 
## Mixed-Effects Model (k = 358; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc 
## -378.4691   756.9381   764.9381   780.4266   765.0524   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.4478 (SE = 0.0350)
## tau (square root of estimated tau^2 value):             0.6692
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   497656476.14
## 
## Test for Residual Heterogeneity:
## QE(df = 355) = 37664401504.9221, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## QM(df = 3) = 271.4606, p-val < .0001
## 
## Model Results:
## 
##                                estimate      se      zval    pval    ci.lb 
## interventiongrazing exclusion    0.1358  0.1250    1.0865  0.2773  -0.1092 
## interventionsoil                -0.7540  0.0472  -15.9653  <.0001  -0.8465 
## interventionvegetation           0.2479  0.0632    3.9228  <.0001   0.1240 
##                                  ci.ub 
## interventiongrazing exclusion   0.3807      
## interventionsoil               -0.6614  *** 
## interventionvegetation          0.3717  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#outcomes
mod.5 <- rma(lrr, var.es, slab= ID, mods=  ~outcome+aridity.index+exp.length -1, data = mdata.all, subset = paradigm == "active")
summary(mod.5)
## 
## Mixed-Effects Model (k = 1102; tau^2 estimator: REML)
## 
##     logLik    deviance         AIC         BIC        AICc 
## -1456.9639   2913.9278   2927.9278   2962.9237   2928.0307   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.7478 (SE = 0.0331)
## tau (square root of estimated tau^2 value):             0.8648
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   105846008972.89
## 
## Test for Residual Heterogeneity:
## QE(df = 1096) = 6497561383370.6289, p-val < .0001
## 
## Test of Moderators (coefficients 1:6):
## QM(df = 6) = 269.7943, p-val < .0001
## 
## Model Results:
## 
##                 estimate      se     zval    pval    ci.lb    ci.ub 
## outcomeanimals    0.1963  0.2040   0.9626  0.3358  -0.2034   0.5961      
## outcomehabitat    0.0193  0.1158   0.1663  0.8679  -0.2077   0.2462      
## outcomeplants     0.6239  0.0789   7.9103  <.0001   0.4693   0.7785  *** 
## outcomesoil      -0.1153  0.1500  -0.7683  0.4423  -0.4093   0.1788      
## aridity.index    -0.0087  0.0042  -2.0854  0.0370  -0.0169  -0.0005    * 
## exp.length        0.0020  0.0003   7.5039  <.0001   0.0015   0.0026  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.6 <- rma(lrr, var.es, slab= ID, mods=  ~outcome+aridity.index+exp.length -1, data = mdata.all, subset = paradigm == "passive")
summary(mod.6) #aridity.index non significant
## 
## Mixed-Effects Model (k = 358; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc 
## -373.2789   746.5579   758.5579   781.7567   758.8006   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.4420 (SE = 0.0347)
## tau (square root of estimated tau^2 value):             0.6648
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   493779723.64
## 
## Test for Residual Heterogeneity:
## QE(df = 353) = 37647543969.8621, p-val < .0001
## 
## Test of Moderators (coefficients 1:5):
## QM(df = 5) = 283.0688, p-val < .0001
## 
## Model Results:
## 
##                 estimate      se     zval    pval    ci.lb   ci.ub 
## outcomehabitat    0.2340  0.1221   1.9161  0.0553  -0.0054  0.4734    . 
## outcomeplants     0.5128  0.1465   3.5002  0.0005   0.2256  0.7999  *** 
## outcomesoil      -0.0521  0.3463  -0.1503  0.8805  -0.7308  0.6267      
## aridity.index     0.0004  0.0064   0.0633  0.9495  -0.0121  0.0129      
## exp.length       -0.0012  0.0006  -1.8689  0.0616  -0.0025  0.0001    . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod.6b <- rma(lrr, var.es, slab= ID, mods=  ~outcome+exp.length -1, data = mdata.all, subset = paradigm == "passive")
summary(mod.6b)
## 
## Mixed-Effects Model (k = 358; tau^2 estimator: REML)
## 
##    logLik   deviance        AIC        BIC       AICc 
## -373.8439   747.6878   757.6878   777.0343   757.8602   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.4407 (SE = 0.0346)
## tau (square root of estimated tau^2 value):             0.6638
## I^2 (residual heterogeneity / unaccounted variability): 100.00%
## H^2 (unaccounted variability / sampling variability):   490927793.15
## 
## Test for Residual Heterogeneity:
## QE(df = 354) = 37647545438.6352, p-val < .0001
## 
## Test of Moderators (coefficients 1:4):
## QM(df = 4) = 283.8777, p-val < .0001
## 
## Model Results:
## 
##                 estimate      se     zval    pval    ci.lb    ci.ub 
## outcomehabitat    0.2397  0.0839   2.8570  0.0043   0.0753   0.4042   ** 
## outcomeplants     0.5185  0.1142   4.5422  <.0001   0.2948   0.7422  *** 
## outcomesoil      -0.0518  0.3458  -0.1498  0.8809  -0.7296   0.6260      
## exp.length       -0.0012  0.0006  -2.0491  0.0405  -0.0023  -0.0001    * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#forest(m1, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))

#forest(m2, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))

#forest(m3, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))

#forest(m4, sortvar = intervention, fontsize = 10, xlim = c(-2, 2))

#t-tests for lrr diff from 0
#mdata %>%
  #split(.$intervention) %>%
  #purrr::map(~sumz(.$lrr)) 
#effect sizes plots
#need better ci estimates

#active plots
ggplot(simple.mdata, aes(intervention, lrr)) +
  ylim(c(-2,2)) +
  geom_point(position = position_dodge(width = 0.5)) + 
  labs(x = "", y = "lrr") +
  coord_flip() +
  geom_errorbar(aes(ymin=lrr-var.es, ymax=lrr+var.es), width=.05, position = position_dodge(width = 0.5)) +
  geom_hline(yintercept = 0, colour="grey", linetype = "longdash")+
   theme(axis.text.x=element_text(face="bold"),
         axis.text.y=element_text(face="bold"),
        axis.title=element_text(size=12,face="bold"),
        strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold"))

ggplot(simple.mdata.2, aes(intervention, lrr, color = outcome)) +
  ylim(c(-2,2)) +
  geom_point(position = position_dodge(width = 0.5)) + 
  labs(x = "", y = "lrr", color = "") +
  coord_flip() +
  geom_errorbar(aes(ymin=lrr-var.es, ymax=lrr+var.es), width=.07, position = position_dodge(width = 0.5)) +
  geom_hline(yintercept = 0, colour="grey", linetype = "longdash")+
  theme(axis.text.x=element_text(face="bold"),
         axis.text.y=element_text(face="bold"),
        axis.title=element_text(size=12,face="bold"),
        strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold"))

#passive plots
ggplot(simple.mdatap, aes(intervention, lrr)) +
  ylim(c(-2,2)) +
  geom_point(position = position_dodge(width = 0.5)) + 
  labs(x = "", y = "lrr") +
  coord_flip() +
  geom_errorbar(aes(ymin=lrr-var.es, ymax=lrr+var.es), width=.08, position = position_dodge(width = 0.5)) +
  geom_hline(yintercept = 0, colour="grey", linetype = "longdash")+
  theme(axis.text.x=element_text(face="bold"),
         axis.text.y=element_text(face="bold"),
        axis.title=element_text(size=12,face="bold"),
        strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold"))

ggplot(simple.mdatap.2, aes(intervention, lrr, color = outcome)) +
  ylim(c(-2,2)) +
  geom_point(position = position_dodge(width = 0.5)) + 
  labs(x = "", y = "lrr", color = "") +
  coord_flip() +
  geom_errorbar(aes(ymin=lrr-var.es, ymax=lrr+var.es), width=.08, position = position_dodge(width = 0.5)) +
  geom_hline(yintercept = 0, colour="grey", linetype = "longdash")+
  theme(axis.text.x=element_text(face="bold"),
         axis.text.y=element_text(face="bold"),
        axis.title=element_text(size=12,face="bold"),
        strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold"))

#ggplot(mdata, aes(lrr, color = intervention)) +
# geom_freqpoly(binwidth = 0.5, size = 2) + 
#  xlim(c(-10, 10)) +
#  labs(x = "lrr", y = "frequency", color = "") +
#  geom_vline(xintercept = 0, colour="grey", linetype = "longdash") +
#  scale_color_brewer()

#ggplot(mdata, aes(lrr, fill = intervention)) +
 # geom_dotplot(binwidth = 1) + 
  #xlim(c(-10, 10)) +
#  labs(x = "lrr", y = "frequency", fill = "") +
#  geom_vline(xintercept = 0, colour="grey", linetype = "longdash") +
#  scale_fill_brewer()

Interpretations

  1. Little to no replication of specific techniques.
  2. Intervention and outcome classification useful and appropriate.
  3. Active versus passive interesting and fascinating.
  4. Use only random effect model statistics because of heterogeneity and conceptual focus of data and synthesis.
  5. Active and passive are not the same - active is net positive and passive is net negative. Amazing finding.
  6. Then subset out each of active passive and run a meta. Found that both for active and passive, vegetation is a path forward for restoration (list techniques in paper) particularly active.
  7. For active interventions aridity index and experiment length are significants in explaining heterogeneity of outcomes. While for passive interventions none of these two moderators were significants.
  8. Soil remediation is powerful active intervention and ignoring passive changes in soil is a dramatic and negative impediment to restoration.
  9. Habitat, animals and soils are challenging to actively (and passively) restore but plants are viable and significant positive outcomes that can be restored.
  10. Passive restoration confirms that plants can recover passively, but soil does not.